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Summary. In the quantitative analysis of Dynamic Contrast-Enhanced Magnetic Resonance Imag- 
ing (DCE-MRI) compartment models allow to describe the uptake of contrast medium with biological 
meaningful kinetic parameters. As simple models often fail to adequately describe the observed uptake 
behavior, more complex compartment models have been proposed. However, the nonlinear regression 
problem arising from more complex compartment models often suffers from parameter redundancy. In 
this paper, we incorporate spatial smoothness on the kinetic parameters of a two tissue compartment 
model by imposing Gaussian Markov random field priors on them. We analyse to what extent this 
spatial regularisation helps to avoid parameter redundancy and to obtain stable parameter estimates. 
Choosing a full Bayesian approach, we obtain posteriors and point estimates running Markov Chain 
Monte Carlo simulations. The proposed approach is evaluated for simulated concentration time curves 
as well as for in vivo data from a breast cancer study. 
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1. Introduction 



Nonlinear regression problems often suffer from parameter redundancy (Seber and Wild 1989). A 



prominent example for nonlinear problems is compartment models, which are used in a variety of 



applications to model the exchange between different compartments ( Schmid 2010 McNally 2008 



Eriksson 1971 Herbst 19631. In this paper, we will concentrate on redundancy issues in complex 



compartment models for Dynamic Contrast-Enhanced Magnetic Resonance Imaging (DCE-MRI). 

DCE-MRI is an imaging technique which allows to image the perfusion in tissue in vivo. After 
injection of a contrast medium, a series of images is obtained. This series shows the uptake dynamics 
of the contrast medium into the tissue over time. For example, in oncology, analysing the dynamics 
of the contrast medium allows to detect tumours, to specify malignancy and type of tumours, and 
to assess the success of cancer therapies (Padhani et al. 2005 Schmid et al. 2009). Typically, 



quantitative analysis of DCE-MRI is based on compartment models which capture the exchange of 
blood (containing contrast medium) between different, well-mixed compartments. With the help 
of differential equations the form of the concentration time curve (CTC) can be analytically de- 
scribed depending on biologically meaningful kinetic parameters, resulting in a nonlinear regression 
problem. 



^Address for correspondence: Volker J. Schmid, Department of Statistics, Ludwig-Maximilians- 
Universitat, Ludwigstrafie 33, 80539 Miinchen, Germany 
E-mail: volker.schmid@lmu.de 



2 Julia C. Sommer and Volker J. Schmid 



1.1. Compartment models 

Compartment models assuming various tissue architectures of different complexities have been 
proposed for quantitative analysis of DCE-MRI data. The simplest and most frequently used 
models are the Tofts model, here also refered to as IComp model, and the "extended Tofts" model 



assuming only the arterial plasma compartment and one interstitial space compartment (Tofts 
1997||Tofts et al.[|1999|). In the Tofts model the observed CTC C t (t) is described by 



C t (t) = C p (t) *^ trans exp(-fc ep t), 



(1) 



with * the convolution operator, C p (t) the arterial input function (AIF), i.e., the concentration 
of contrast agent in the blood plasma, K trans the transfer rate from blood plasma to extracellular 
extravascular space (EES), and k cp the rate constant for transfer between EES and plasma. 

Tumour tissue is often heterogeneous (Schmid 20101 and the Tofts and extended Tofts models 
fail in describing its observed uptake dynamics (Schmid et al. 2009). Therefore, several authors 



propose more complex models to describe perfusion in tissue. For example, the two compartment ex- 



change model (2CXM) has seperate compartments for arterial plasma and interstitial plasma ( Brix 
et al. 20091 Sourbron and Buckley 2011). Multi-compartment models allow for two to three ki 



netically distinct tissue compartments to describe CTCs on a region of interest level (Port et al. 



1999). Here, we use two tissue compartments (2Comp model) with the same architecture as Port 



et al. ( 1999 ) to model CTCs on a per voxel level: 



C t (t) = C p (t) * (Af ans exp(-fc cpi t) + ans exp(-fc ep2 t)) 



(2) 



With the 2Comp model on a voxel level CTCs in heterogeneous tissue can be more adequately 



described, especially at tumour margins (Karcher and Schmid 2010) and tissue heterogeneity be- 



comes accessible. This is important because tissue heterogeneity is diagnostically informative. The 
2Comp model is introduced in more detail in Section |2.1| and compared to other compartment 
models in Section |2~21 



1.2. Parameter redundancy 

Though necessary, increased model complexity brings challenges. With increased complexity model 
parameters can become redundant and in this case cannot be stably estimated. 

Parameter redundancy (or non-identifiability) is frequently encountered in nonlinear regression 
problems. In contrast to identifiability problems occurring in standard linear regression models, 
problems occurring in nonlinear regression models are often such that they cannot be eliminated 



by optimal design (Gilmour and Trinca, 2012), reparametrisation, nor reduced error variance. This 



is the case when different nonlinear functions corresponding to different models or parameter con- 
stellations are too similar. Seber and Wild describe such an example for a sum of two or three 
exponentials ( Seber and Wild 1989 p. 119): 



Even though these functions are so different, the curves are visually indistinguishable. 
This approximate lack of identifiability is called parameter redundancy by Reich [1981]. 
Models such as this give rise to bad ill- conditioning, no matter where the x's are placed — 
a common problem with linear combinations of exponentials. 



We will discuss the redundancy problem for the 2Comp model in more detail in Section 2.3 



As a solution to the redundancy issues we suggest to regularise the parameter space by spatially 
smoothing the parameter maps. We propose a Bayesian framework in order to determine the 
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kinetic parameters of the 2Comp model on a voxel level, using prior biological knowledge on the 
parameters and accounting for the spatial structure of the image. The approach makes use of the 
intrinsic spatial information given by the voxel structure of the image. Spatial regularisation is 
done by Gaussian Markov random fields (GMRF) as priors on the kinetic parameters as previously 



proposed by Schmid et al. (2006) and Kelm et al. (2009) for the IComp model 



1.3. Outline 

In this paper we aim to analyse to what extent the spatial regularisation helps to obtain stable 
parameter estimates in the 2Comp model. We choose a full Bayesian approach and obtain poste- 
riors and parameter estimates by running Markov Chain Monte Carlo (MCMC) simulations. The 
proposed model is evaluated for simulated CTCs and data from a breast cancer study. 

An advantage of the Bayesian framework is that the posterior can still be computed in the 
case of parameter redundancy; however, the redundant parameters will have multimodal marginal 
posteriors. We find that assuming spatial smoothness on the exponential rates is an efficient way 
to regularise the parameter space and to make parameters of a 2Comp model identifiable. As a 
result we obtain parameter estimates at a voxel level that are more stable. Hence, one can describe 
heterogeneity of the tissue without losing spatial information. 



2. Model 



The observed contrast concentration Yij at time tj, j = 1, T in voxel i = 1, N can be described 
by the theoretical concentration time curve Ct (t) depending on the voxel-specific kinetic parameters 
(jf plus Gaussian noise, i.e., 



Y, 



N {CtW-JiWi 



(3) 



with erf the voxel-specific variance of the Gaussian noise (Schmid et al. 2006). The form of Ct and 
the number of kinetic parameters is fixed, given a specific compartment model. 



2. 1. Two tissue compartment model (2Comp) 











t 2 







plasma 

Fig. 1. Two tissue compartment model 

Even at a voxel level, tissue can be heterogeneous. That is, there may be two tissue compart- 
ments with different kinetic properties that exchange with plasma at distinct rates k ePl and k ep2 
(see Fig. [I]). Hence, we assume a model with two tissue compartments (2Comp). We label the 
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compartments such that fc cPi < k cp2 . This means that the first compartment is the slow and the 
second compartment is the fast exchanging compartment. 



As k 



Kl r ^ ns /v tk (Tofts et al. 1999), the volume v tk of tissue tk per unit volume of tissue 



can be computed as v tk = Kl rans / k ep . The changes in tissue concentrations are given by 



v t ±C t ,{t) = Kf 
d - t C t2 (t) =Ktf™ 



[C p {t)-C tl {t)] 
«t2 3t°*2W = «2'""' [c P (t) ~ CtS)\ 

The solution of these differential equations is given by 



C tk (t) = C p (t) 



j^trans 

* — exp 

v t k 



Kt 



-t 



(4) 



(5) 



for k=l,2 denoting the different tissue compartments. The total (observable) concentration is then 
given as C t = v tl C tl + v t2 C t2 by 



Ct(t) 



E 



C p (t)*Kl Ians exp(-k ePk t). 



(6) 



The AIF describes the input of contrast agent through the blood stream. As suggested by Tofts 
and Kermode ( 1991 ), we use a bi-exponential function of the form 



C p (t) = Dy^aiexp(-mit) 



(7) 



with dose D and values a\ = 3.99 kg/l, a2 = 4.78 kg/l, mi — 0.144 min l , m^ = 0.0111 min . 

We use an exponentia l paramctrisation that insures the rate and transfer constants to be positive 
(see Schmid et al. 2006 and references therein): 9 l k = log (kl Pk ^J , l\ = log ^* rans ' 1 ^ ] for = 1, 2. 



2.2. Relation to other compartment models 

In the proposed 2Comp model, the observed concentration Ct(t) is described by an impulse response 
function (sum of two exponentials) convolved with the AIF, see equation ([6]). In the 2CXM the 
interstitial space and the interstitial plasma are modeled with seperate compartments. Though 
explained by different compartment designs, the impulse response of the 2CXM is also a sum of 
two exponentials. Hence, the 2Comp model and the 2CXM lead to the same nonlinear regression 
problem. We prefer to use the 2Comp model due to the charming fact that the impulse response 
is directly expressed by interpretable parameters i^J; rans , K^ aris , k ePl and k cp2 . In contrast, in 
the 2CXM the impulse response is expressed by auxiliary variables (called F + , K + and KJ) 
which are complicated functions of interpretable quantities — see Lemma 3 of |Sourbron and Buckley] 

pm] ). 

For the case that the exchange rates are the same, fc epi = fc ep2 , or when one of the tissue volumes 
vanishes, Vt t = or v t2 = 0, the impulse response reduces to a single exponential and the 2Comp 
model corresponds to the standard Tofts model, also referred to as IComp model here. 

For the case that one exchange rate becomes infinite, k ep2 — oo, the observed concentration is 
of the form 

C t (t) = v t2 C p (t) + C p (t) * Kl ians exp(-k cpi t). (8) 

In this case the second tissue compartment takes the role of an interstitial plasma compartment, 
Ct, = C p , and the 2Comp model corresponds to the extended Tofts model. 
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2.3. Redundancy issues in the independent voxelwise model 

In a voxelwise approach, the CTCs of all voxels are fitted independently from each other. Similar 



to the voxelwise Bayesian IComp model evaluated in Schmid et al. (20061, we impose Gaussian 
priors on the logarithmic rate constants 61 

and on the logarithmic transfer constants 7^, 

independently for all voxels i = 1, . . . , N with fixed precisions Tg k — T Jk — 1 and fj,$ 1 — = /i 7 , = 
0, fig 2 = log(5) . With this prior, all rate and transfer constants kl Pk and ^f* rans,z remain positive. 
Rate and transfer constants of the first compartment with a priori probability of 99.86% do not 
exceed 20 min . The dynamics in the second compartment is assumed to be faster with a priori 
expected k\ p ^ values of five. 




Fig. 2. Similar CTCs for two different parameter vectors. Black line: CTC described by fc ePl = 2.07, fc eP2 = 
2.07, K x { ans = 0.55, A"2 ans = 0.15 (can as well be described by only one compartment with k ePl = 2.07 and 
^trans = Q 7 y Grey || ne: CTC described by fcepi = 2 .i9, fc e p 2 = 5.02, iff 3 " 5 = 0.62, _R-!, rans = 0.24. The dashed 
lines show the contribution of the first compartment and the dotted lines those of the second compartment. 



This independent voxelwise model leads to unstable estimates due to redundancy issues ( Karcher 



and Schmid 2010). Obviously, redundancy issues occur when the contribution of one compartment 
vanishes. However, they may as well occur when the exponential rates are too similar. There 
are theoretical results on parameter redundancy in sum of exponentials models (Seber and Wild 



1989), however, a generalisation for the case of convolved exponentials is tricky. In Reich (1981) 



a redundancy measure was used to show that parameters in a sum of two exponentials model are 
highly redundant if the exponential rates differ by less than a factor of five. Even though this result 
does not directly transfer to the convolved exponentials given in equation ([6]), this results still helps 
to understand parameter redundancy in the 2Comp model. In Fig. [2] we show an example for data 
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simulated from a two-compartment model where the exponential rates differ by a factor of four. 
In this case, a solution from a compartment model with only one tissue compartment exists that 
fits the observed concentration reasonably well. Therefore, several quite distinct parameter vectors 
describe very similar CTCs, meaning that parameters are redundant. 



2.4. Spatial regularisation 

In order to solve the problem of redundancy, we propose to use spatial regularisation. In the 
following, we will introduce a spatial prior which accounts for the spatial information intrinsic in 
the DCE-MR images. 

In the proposed spatial model, we assume that rate constants fc* Pfc vary smoothly in space and 
hence that the exchange properties of each tissue compartment are rather smooth. In contrast, the 
contribution of differently behaving compartments in each voxel is assumed to be quite flexible, 
meaning that the tissue volumes v\ k — ^* rans ' 1 /fc* Pfc may vary unsmoothly from voxel to voxel. 

Then the transfer rate as a product of rate and volume j^^ ans ' 1 = k^ v\ k inherits the spatial 
smoothness of k l epk , but is less smooth due to varying v\ k values. The spatial smoothness of the 
kinetic parameters is modelled using a Gaussian Markov random field on its logarithmic transforms 



61, 7 J. (Rue and Held 2005 Schmid et al. 2006) 



We use a neighbourhood structure where adjacent voxels are neighbours, that is, each voxel 
has four neighbours unless it lies at the edge of the image. From this, a prior distribution can be 
defined by assuming a Gaussian distribution on the differences of neighbouring logarithmic rate 
and transfer constants: 

61 - el\re h ^ N(0, (re,)- 1 ) for ir^j 

and 

7* --^K^JVCO,^)- 1 ) far 

This spatial prior on the logarithmic rates leads to smooth parameter maps of kl Pi and fcg Pa , 
j^trans,i an j j^trans^ However, a priori we expect much smoother maps for fc* pi and fc* p2 and 
less smooth maps for ^ rans ^ anc j X 2 raJls ' 1 . jj ence; we use Gamma priors on the precisions Tg k ~ 
Ga(ag k ,bg k ) for fc = 1,2 with ag 1 = ag 2 = 1000 and bg 1 — bg 2 — 1 and T lk ~ Ga(a 7k ,b 7k ) for 
k = 1, 2 with a 7l = a 72 = 0.0001 and 6 7l = 6 72 = 0.01. 

Furthermore, we assume the noise variance to be the same of = a 2 for all voxels i. For the 
observation variance we assume an Inverse Gamma prior af ~ IG(a, b) with a and b such that the 
a priori expected SNR corresponds to values typically observed in breast tumour DCE-MRI data 
(SNR typically ranges from 10 to 20). We choose this prior to be more informative with increasing 
number of voxels. 



2.5. Implementation 

We implemented the proposed spatial 2Comp model extending the R-package "dcemriS4" ( Whitchei 



and Schmid 2009 20111. For each voxel, we simulate from the posterior of the model param- 
eters with a MCMC algorithm (Gilks et al. 1996). Starting with random values, the voxels 



are subsequently updated in random order. More precisely, starting values are drawn from uni- 



form distributions per voxel with 



C/[0,1], 



= 1-vf 



restart 
K ep 2 



J7[1.75, 5.25]. The logarithmic rate and transfer constants 6\, 9^, 7i, 



and kf u art 

e Pl 



U[0. 1,0.3], 
72 are updated with 
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Metropolis Hastings steps with random walk proposal. Gibbs update steps are used for the hyper- 
parameters l/<7 2 , Tg k and T* fc , as its full conditionals are Gamma distributions that can be sampled 
from directly, see Appendix |A.1| 

The proposal variances of the random walk proposals are tuned such that the Metropolis- 
Hastings acceptance rates are approximately 20%. After a burn-in of 5,000 iterations, 5,000 iter- 
ations are performed with every third sample saved. For parameter point estimation we use the 
median of the MCMC sample. 



2.6. Measure of model complexity 

We suggest to use the number of effective parameters pn as a measure of model complexity and, 
hence, tissue heterogeneity. The number of effective parameters po is calculated as the differ- 
ence of the posterior median of the deviance and the deviance evaluated at the posterior median 
value ( Spiegelhalter et al. 



2002). It was introduced for the calculation of deviance information 



criterion (DIC) and is typically used to asses model fit and complexity in Bayesian hierarchical 
models. The DIC is defined as the posterior median deviance plus the number of effective parame- 
ters po (ISpiegelhalter et al. 1 120021). 



However, as 



Spiegelhalter et al. (20021 point out, po can become negative in cases where the 

This is certainly the case when dealing 



posterior mean or posterior median is a poor estimator 
with multimodal posteriors due to parameter redundancy. Hence, we can detect redundancy issues 
by looking at the pn- 

Here, although we are dealing with a joint model for all voxels, i.e., a model for the whole 
image, we will also compute a voxelwise pd using the deviance in each voxel. The deviance is 
D((j) l ,a 2 ) = — 2 * Z(0*,cr 2 ) where l{(j) l ,(j 2 ) is the log-likelihood function given in Appendix A.l 
is evaluated at the posterior median values of 



It 

and a 2 in order to calculate the deviance of the 



median and it is evaluated at each sample value of (f> 1 and a in order to calculate the median 
deviance. This allows to assess the model complexity per voxel and, hence, the tissue heterogeneity. 



3. Simulation study 

3. 1 . Simulation setup 

In order to evaluate the proposed voxelwise and spatial 2Comp models, we simulated a DCE-MR 
image of 25 x 25 voxels with different parameter combinations in a two tissue compartment model. 
The parameter configuration is given in Fig. [3] For three blocks of different size, we simulated CTCs 
from a true 2Comp model, i.e., a mixture of two tissue compartments with very different exchange 
rates fc opi = 0.2 and k cp2 = 4. Both compartments contribute equally with volumes v tl = v t2 = 0.5. 
In the lower left block, the exchange rate k epi varied smoothly from 0.2 in the middle to 0.5 at the 
corner. 

For two blocks, we simulated CTCs from a IComp model. One of those blocks is described by 
a tissue compartment exchanging rather slowly with plasma at rate fc ePi = 0.2 and having volume 
Vt 1 — 1 (the fast exchanging compartment with k cp2 = 4 has no contribution, i.e. Vt 2 — 0). For 
the other block, the exchange with plasma is rather fast: k ep — 4, vt 2 — 1 (the slow exchanging 
compartment k ePi = 0.2 has no contribution, i.e. Vt ± = 0). 

Within each block, uniformly distributed noise U[0.8, 1.2] was multiplied to the parameters A: ePl , 
k cp2 , if* rans and i^™ 13 P er voxel and the corresponding CTC was computed from these kinetic 
parameters. Gaussian noise was added to the simulated CTCs with standard deviation a — 0.05. 
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kep t = 0.2 
kep 2 = 4 



(v„=0) 



v„ =0.5 
v„ = 0.5 



IComp 



IComp 



IComp 



IComp 



kept = 0.2 



fep, = 0.3 



kepi = 0.4 



kep 2 = 4 



(v„ = 0) 



kep x = 0.5 



Fig. 3. Sketch of simulation design 



3.2. Results 

As discussed above, when fitting a model with two tissue compartments one often deals with iden- 
tifiability issues. In these cases, the model is overparametrised and one observes unstable estimates 
of all parameters. The Bayesian approach allows to evaluate the posterior anyway; however, in the 
case of redundancy, the marginal posteriors typically are multimodal. 

For example, in the blocks simulated from a IComp model parameter estimates are redundant 
and hence unstable in the lower right of the simulated image. For one of these voxels, Fig. [4] 
depicts the marginal posteriors of volume fractions v tk and rate constants fc ePfc . With the voxelwise 
approach, the posteriors are multimodal and hence there is no good point estimator. In comparison, 
the spatial approach produces unimodal posteriors and good point estimates can be gained by 
computing the median of the MCMC sample. In contrast to the voxelwise model, the contribution 
of the two compartments are well separated (fc oPl and k cp2 samples are not too similar) and the 
estimated volume of the first compartment is close to zero. 

Using the point estimates for the kinetic parameters we obtain estimated CTCs and refer to 
them as fit. In Fig. [5] we compare the fit of the IComp model, the voxelwise, and the spatial 2Comp 
model. For a curve simulated from a IComp model (Fig. [5] (a)), the fit of the spatial 2Comp model 
is similar to the fit of the IComp model. However, the voxelwise 2Comp model fails to adequately 
fit the curve due to redundancy issues. The SSE is about 0.14 for the voxelwise 2Comp model and 
about 0.1 for the IComp model as well as the spatial 2Comp model. For a curve simulated from 
a 2Comp model (Fig. [5] (b)), both the spatial and voxelwise 2Comp models clearly outperform the 
IComp model with similar fits. Here, the SSE is about 0.12 for the spatial and voxelwise 2Comp 
models and 0.2 for the IComp model. 

In Pig.|6jthe sum of squared errors (SSE) per voxel are depicted for the IComp and the voxelwise 
and the spatial 2Comp model. Considerable differences in SSE for the IComp model compared to 
both 2Comp models can be observed for the three blocks simulated from a true 2Comp model. The 
voxelwise and spatial 2Comp models have similar SSE, with increased SSE in the voxelwise model 
for voxels with multimodal posteriors. These differences cannot be distinguished at this scale and 
were shown for a specific curve above (Fig. [5] (a)). 

For voxels with multimodal posteriors, the estimates of pjj are not meaningful. In the voxelwise 
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Fig. 4. Posterior marginal densities for curve simulated from 1Comp. Left: Voxelwise 2Comp fitted. Right: 
Spatial 2Comp fitted 



2Comp model estimated pjj values are often negative due to parameter redundancy. In contrast, 
Pd values are always positive in the spatial 2Comp model as the posteriors are unimodal. Values 
range between 0.5 and 1 with median 0.68 in the IComp blocks. In the 2Comp blocks, pn values 
between 0.8 and 1.6 with median 1.2 show increased tissue heterogeneity. 

In Fig. [7]we show the parameter maps for the estimates of £; 0Pl , fc 0P2 , if* rans and Kl ians from the 
voxelwise and the spatial model as well as the true underlying parameter values. As the voxelwise 
approach leads to unstable estimates, the estimation results differ strongly from the true underlying 
values. Especially for the voxels simulated from a IComp model, the voxelwise 2Comp model leads 
to unstable estimates. For instance, for voxels in the lower right simulated from a true IComp 
model — with vt 1 = — k cPi is overestimated, Kf ans is overestimated, and consequently K^ &ns is 
underestimated. Compared to the voxelwise model, the spatial model leads to smooth parameter 
maps that can be interpreted more intuitively and to stable estimates that better match the true 
underlying parameter values. 

In the spatial model, the parameter maps for fc GPi and fc ep2 are smooth and the estimates match 
the true underlying values quite well. There is some oversmoothing such that the higher fc GPi values 
in the lower left corner are underestimated and as a consequence also the corresponding fc cp2 are 
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Fig. 5. Simulated data and curve fits for 1Comp model, spatial and voxelwise 2Comp models, (a) Data 
simulated from a 1Comp model (parameter constellation as in the lower right corner), (b) Data simulated 
from a 2Comp model (parameter constellation as in the upper right corner). 



underestimated. The estimates of K\ rans and K^f 3 * 15 perfectly match the true underlying values. 
For the blocks simulated from a IComp model either the Kl rans estimate or the J^* rans estimate 
becomes zero. Like this, model redundancy is avoided and the posteriors are unimodal. 



4. Breast cancer study 

4. 1 . Data description 

To evaluate the clinical use of our approach we use a previously analysed DCE-MRI study on breast 



cancer (Schmid et al. 2006). The dataset consists of twelve patients with breast tumour. The scans 
were acquired once before and once after six cycles of chemotherapy. The drug is expected to stop 
the process of angiogenesis, i.e., to lower the elevated blood flow to the tumour, hence, to lower 
the kinetic parameters K tYans and fc op . As clinical evaluation, tumours were removed at the end of 
therapy and the response to therapy was evaluated histologically. Six of the twelve patients were 
identified as responders, the other six as nonresponders. Informed consent was obtained from all 
patients and the study was acquired in accordance with recommendations given by |Leach et al. 



(2005) 



The scans were acquired with a 1.5 T Siemens MAGNETOM Symphony scanner, TR = 11 ms 
and TE = 4.7 ms. Each scan consists of three slices of 230 x 256 voxels, but only the central slice 
was used in our analysis. A dose of D = 0.1 mmol/kg body weight Gd-DTPA was injected at 
the start of the fifth acquisition using a power injector. Regions of interest cover the tumour and 
surrounding normal tissue. 



4.2. Results 

In Fig. [8]the parameter maps from the spatial 2Comp model are shown for pre- and post-treatment 
scans of patient 4 (nonresponder to therapy) and patient 6 (responder). Similar to the results 
of the simulation study and in accordance with the prior assumptions, the estimated parameter 
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Fig. 6. Evaluation of model fit: sum of squared error (SSE) and p D 



maps for the DCE-MR images are quite smooth for the exponential rates fc oPl and fc e p 2 , whereas 
the K\ lans and ifJj rans estimates show more spatial variation. The contribution of the second 
compartment vanishes [K^ 3 '™ close to zero) in healthy tissue. In those regions, the IComp model 
suffices to describe the observed uptake dynamics, meaning that the tissue is homogeneous there. 
Interestingly, tissue inside of the tumour is often homogeneous as well. The second compartment 
has nonzero contribution and improves the fit of observed CTCs at tumour margins and in parts of 
the surrounding tissue. In those regions, the tissue is heterogeneous as both the slow and the fast 
exchanging compartments contribute to the uptake dynamics. Larger pp values and improved fit 
compared to the IComp model reflect this heterogeneity. 

For patient 4 the parameter maps for the pre-treatment scan depict high fc opi and k cp2 values 
as well as high if' rans and K2 lans values for a large tissue region. Post-treatment, the kinetic 
parameters fc ePl , fc op2 , K\ lans and i^| rans have higher values, but the tissue region with increased 
blood flow becomes smaller and more dense. Reduced tumour volume could easily be misinterpreted 
as treatment success. Here, the pu map contains additional information that might help to assess 
treatment success. For patient 4 voxelwise po values are even higher in the post-treatment scan. 

For patient 6 parameter maps of fc ePl and fc eP2 are quite smooth. The contribution of the 
second compartment — — is close to zero inside of the tumour and in surrounding healthy 
tissue, and it is not vanishing at tumour margins and in surrounding tissue. After treatment, 
the number of voxels where the second compartment contributes decreases notably. For patient 6 
tumour margins and extensions around the tumour are heterogeneous and better described with 
the aid of an additional second compartment. Both the pu values and the size of the tissue region 
with increased pu decrease after treatment. 



1 2 Julia C. Sommer and Volker J. Schmid 
5. Conclusions 



In this paper we have discussed redundancy issues of a specific nonlinear regression problem, namely 
the estimation of kinetic parameters in a two tissue compartment model for DCE-MR images. With 
a spatial prior we have regularised the parameter space and made the parameters identifiable. With 
this prior, a 2Comp model can be fitted at a voxel level and CTCs in heterogeneous tissue, especially 
at tumour margins, can be described better than with the standard IComp model. For CTCs that 
are adequately described by the IComp model, the estimates of one of the compartment volumes 
is close to zero. Like this and in contrast to a voxelwise approach parameter estimates are stable 
and easy to interpret. 

Confronted with redundancy issues, modelling with compartments requires trade-off between 
too simplistic models and overparametrised, redundant models. Complexity is determined by the 
number of compartments as well as the level of spatial resolution. Spatial regularisation offers 
a solution that can be applied in other fields as well, for instance in the quantitative analysis of 
positron emission tomography (PET) and single-photon emission computed tomography (SPECT) 
images. In PET and SPECT neuroreceptor imaging studies the kinetics of ligand uptake in the brain 



is described with the help of compartment models (Slifstein and Laruelle 2001). When estimating 



receptor parameters one copes with similar identifiability issues encountered in DCE-MRI analysis. 

We have proposed po as a measure that contains additional information about the heterogeneity 
of the tissue whereas the kinetic parameters contain information about the uptake dynamics only. 
We find it interesting that estimates of the effective number of parameters pjj rarely exceed values 
of 1.5, even for CTCs simulated from a 2Comp model with four kinetic parameters. 

In summary, we proposed and evaluated spatial regularisation for two-compartment models, 
allowing a more comprehensive insight into tissue perfusion, in particular in heterogeneous tissue. 
Spatial regularisation allows to overcome the redundancy issues by "borrowing strength" across the 
tissue of interest, and hence allows to fit complex compartment models even on voxel level with low 
signal-to-noise ratio. Additional clinical studies should be performed to further explore the clinical 
potential of this model. 
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A. APPENDIX 

A. 1. Likelihood and full conditionals 

The log-likelihood depends on the voxel-specific kinetic parameters <p l and the inverse noise variance 
~ _ i . 



T 

K^T,) = |log(2 7 rr e )-ir £ ^(r iJ -C t (^,i J )) 2 . (9) 



.7=1 



In the spatial model, the full conditional distribution of the logarithmic rate constant in voxel 
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i, 9 l k , given the logarithmic rate constants of all other voxels, 9 k % , 

p(0j|V,TOaexp(-^ E (fk-Oi) 2 ) 

jed(i) 

depends only on those of its direct neighbours for k — 1,2. Here, d(i) denotes the set of direct 
neighbours of voxel i. The full conditionals of the logarithmic transfer constants jI have the same 
form. 

Let eij = Yi.j — Ct{<fi l ,tj) denote the random noise terms. Then, the full conditional of t c is 
r e \- <~ Ga(a + ^l£-,b + \ J2iLi Sj=i e 1j) f° r the spatial model. The full conditional of the precision 

< ^ < IV ~ Ga(a g + b e + \ £ . £a(i) (rj i - rtf . Similarly for rj 2 , and r^. 
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Fig. 7. Parameter maps for simulation study: Voxelwise (left column) and spatial fit (middle) of 2Comp model 
for true underlying values (right) 



Spatial two tissue compartment model for DCE-MRI 



patient 4 patient 6 

pre post pre post 




Fig. 8. Spatial 2Comp model: parameter maps for patients 4 and 6 pre- and post-treatment 



